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We propose new approach to numerical study of quantum spin systems. Our 
method is based on a fact that one can use any set of states for the path integral as 
long as it is complete. We apply our method to one- dimensional quantum spin sys- 
tem with next-to-nearest neighbor interactions. We found remarkable improvement 
in negative sign problem. 
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Section 1 Introduction 

Recently quantum spin systems have obtained much interests among people in 
various fields. One reason for it is that Haldane found characteristic property of the 
quantum spin system which is difficult to imagine in the classical case [|I|. Another 
reason comes from the possible relationship between the antiferromagnetic system 
on a 2-dimensional square lattice and the high materials [0. 

One powerful tool to numerically investigate quantum spin systems is Monte 
Carlo approach using the Suzuki- Trotter formula 0. Study through this method 
has brought us very intriguing results on the ferromagnetic system. If one applies 
this method to the antiferromagnetic system, however, one often encounters the 
so-called negative sign (N.S.) problem. N.S. problem, which becomes more serious 
on larger lattices, makes it very difficult to get statistically meaningful results in 
numerical calculations. 

Let us explain what the N.S. problem is. To this purpose we show how the Monte 
Carlo method is applied to the quantum spin 1/2 system whose Hamiltonian is 

where cTj = {af^a^^a^) represents Pauli matrix on site i of a lattice and sum runs 
over all nearest neighbors. The partition function Z is calculated by 

Using the Suzuki- Trotter formula it can be written by the classical partition function 
as follows. 

Z = lim tr|(e-^^i/"e-^^^/")"} 

n 

= i™o ^ n < I e"^^^/" I a[ >< a[ \ e'^^^/" | a^+i >, 

H = Hi + H2. 

Here | ctj > represents a state of the system, {| ctj >} a complete set of states 

1 = X! I X I' 
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and H is divided into two parts so that every term in Hi {H2) commutes with other 
terms within that partial Hamiltonian. Since the expectation value of e~^^^'^ is 
a real number we can apply Monte Carlo methods to this partition function if n, 
which is called trotter number, is kept finite. 

There would have been no N.S. problem if positivity of the expectation value were 
always guaranteed. However in several cases of much physical interest, antiferro- 
magnetic quantum spin system for example, we do have some negative expectation 
values. If total product of expectation values over one configuration can become 
negative, then the N.S. problem may occur. In order to numerically calculate some 
physical quantity (A) one should subtract contributions of negatively signed config- 
urations, A_, from those of positively signed ones, A^, namely. 




where Z+(Z_) is number of configurations with positive (negative) weight. The 
result would suffer from serious cancellation when Z_ ~ Z+. 

Note that negativity of expectation values does not always bring the N.S. problem. 
On square lattices one does not need to worry about it even in the antiferromagnetic 
case because the total product of expectation values is always positive. On the 
triangular lattice, on the contrary, the N.S. problem is quite serious. 

In this paper we describe our prescription to obtain meaningful numerical results 
for systems to which the conventional formulation is useless because of the N.S. 
problem. Our strategy is to find rearranged basis of the system with which no serious 
cancellations take place. We do not mean that the N.S. problem is completely solved 
by our approach, but we show our method is successful to improve the numerical 
results when applied to a quantum spin system. 

In following sections we concentrate ourselves to one dimensional quantum spin 
1/2 system with the next-to- nearest neighbor interactions, which is the simplest 
one among systems suffering from serious N.S. problem. We show what choice 
of complete set of states is more suitable in this case to carry out Monte Carlo 
simulations. In section 2 we present our method. Numerical results will be given in 
section 3 and final section is devoted to discussions on technical problems. Physical 
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properties of the system will be discussed in other paper. 

Section 2 Method 

Hamiltonian of the quantum spin 1/2 system with next-to-nearest neighbor in- 
teractions on a one- dimensional chain is 

J ^ ^ ^ 

^ i=l 

where N is number of sites on the chain and (Tat+i = oi (periodic boundary con- 
dition). The system is ferromagnetic (antiferromagnetic) for positive (negative) J. 
We imply the coupling constant between next-to-nearest neighbors is equal to that 
between nearest neighbors. It is partly because for simplicity and partly because 
this case brings most serious N.S. problem for negative J. The partition function Z 
is given by Z = tr{e~^^). 

First let us describe the conventional approach. State on each site is represented 
by 2;-component of the spin, namely up and down, or -|- and — . In this representation 
states of the system are given by 

I a >=| si, 82,83, ...,8n >, 

where Si = + or — . The identity operator is then 

i = I Sl, S2, ■-, sn >< Si, 82, 8n I ■ 

{Si} 

To use this identity operator in the Suzuki- Trotter formula, we divide the Hamilto- 
nian into four parts as schematically shown in Fig. 1(a). 

n— »oo 

J ^^"^ 
^ i=l 
J 

H2 = -— X] ^2i^2i+l) 
^ i=l 

J ^''^ 

H3 — —— X](^4i-3^4i-l + '?4i-2'?4i), 
i=l 
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J ^-^^ 

H4 = —— y^(^4i-1^4i+l + <?4i^4j+2)- 

With these complete sets and partial Hamiltonians we obtain the partition function 
Z^^ which we use in Monte Carlo calculations, 

n 

= E n < I e"^^^/" I X I e-^^^/" I a'; > 

{aj,a'.,a'!,a'!'}j=l 

x<a]\ e-/^^^/" I >< I e-'^^^/" | >, 

where suffix j numbers sites along the trotter axis and a^+i = ai. 

Now we turn to our method. First we will show our idea using spin system of two 
sites which we denote site a and site b. The Hamiltonian of this system is simply 
Hab = —Joa<^b- In the conventional approach to this two-spin system the complete 

set is given by {| +„, +b >, | +a, —b >, \ —a, +b >, I —a, —b >}• Instead of this set we 
can adopt another complete set composed of diagonalized states of Hab, 

{| 1 >,| © >,| e >,| -1 >}, 

where 

I 1 > = | +a,+b >, 

1 

I ® >= ^(1 +a, -b>+{\ -a, +b >), 

I e >= -^(1 +a, -6 > -(I -a, +b >), 
I -1 > = | -a,-b > ■ 

The partition function is then easily calculated using eigenvalues Ei, E^, Eq and 
E-i, 

Zab = trie-f^^-^o) 

We apply this diagonalization to the spin system of N sites. For this purpose we 
rewrite the Hamiltonian in the following form 

H — — — y^X^a,i^a4+l + ^b,i^b,i-^\ + ^a,i^b,i + ^b,i^a,i+l), 



5 



where we denote odd and even sites with suffix a and b, respectively, 



We employ the complete set where the operators (^a^iCTb^i are diagonalizcd. Then 
each state of the system is represented by combinations of N/ 2 two-site diagonalizcd 
states, 

I Q; >=| Si, S2, <S'jv/2 >, 

where stands for diagonalizcd state on the z-th spin pair, namely Si = Ij, ©j, 
Qi or — Ij. Then it becomes necessary to divide the Hamiltonian into "odd" and 
"even" parts 

H — Ho + He, 

where, as shown in Fig. 1(b), 

J 

Ho — —— ^(o^a,2i-l0^a,2i + (?fe,2i-l'?6,2i + (?b,2i-l0^o,2i 
^ i=l 

+ -Ca,2i-lC6,2i-l + -0'a,2iO'b,2i) , 

J 

He = —-^ X/(^«.2i'^a,2i+l + '?6,2i^6,2i+l + <?6,2i^a,2i+l 
^ i=l 

U ^ 1- - ^ 

+ -^Cra,2i-lO'b,2i-l + -Ca,2iC6,2ij- 

The partition function with these partial hamiltonians is 

Z= lim ir{(e-^-^°/"e-^-^^/")"|. 

n— »oo 

Inserting identity operators made of a's between the exponents we obtain partition 
function Z^'' to use in our numerical study, 

n 

Zm- E n < I e"^^"^" I X I e-^^^/" I > . 

{aj,a'.}j=l 
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Section 3 Numerical Results 

In this section we present Monte Carlo results we obtained using the new partition 
function stated in section 2. In order to figure how much quantitative improvement 
we can make, we also show some results from the conventional approach. 

Monte Carlo simulations are carried out on a two dimensional lattice, one of 
which is space direction of size N and another is the trotter direction having 2n 
(4n) sites in our (conventional) method with the trotter number n. 

In each simulation with new partition function typically one thousand con- 
figurations are generated in order to reach thermal equilibrium and ten thousand 
configurations after thermalization are used to measure system's energy, E, and 
ratio of negatively weighted configurations to total configurations, P, 
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Here one new configuration is obtained by updating spins on all sites using the heat- 
bath method. In each update we make simultaneous spin fiips on sites along the 
smallest closed path (local fiips) and along the trotter axis (n-direction global fiips) 
according to their probabilities. From technical reasons we did not include local and 
global fiips which change 2;-component of Si of i-th spin pair by 2 . We also made 
global fiips along the space axis (x-direction global fiips) available, but we excluded 
them in simulations presented in this paper because we have found their effects on 
system's energy are negligible. 

In simulations with the conventional partition function we employ Metropolis 
algorithm. Here, in addition to the local flips and usual n-direction global flips, we 
found it necessary to make a kind of twisted n-direction global flips shown in Fig. 
1(c) so that enough phase space is guaranteed. We carry out twenty thousand 
updates for the thermalization, one hundred thousand for the measurement. 

Let us flrst present results for the ferromagnetic system. In this case the con- 
ventional method should be effective since no N.S. problem exists. Fig. 2 shows 
results on an = 8 chain by the conventional method as a function of temperature 
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T, which is defined by 1/ J (3, together with exact values obtained by diagonahza- 
tion methods. We see conventional method works very well. Our method is also 
successful, but we have no special reason to prefer it [Q. 

Next we will show results for the antiferromagnetic system. In Fig. 3(a) we 
present ratio P versus T with A^=8 chain by both methods. As we see in the figure, 
even when n = 2 the ratio P with the conventional partition function increases quite 
rapidly as T decreases. It is therefore difficult to obtain statistically meaningful 
results by this method. By our new method, on the other hand, P remains much 
smaller for larger n and for lower T so that the N.S. problem is less serious. Fig. 
3(b) plots system's energy on = 8 and n = 2, 3 and 4 lattices in comparison with 
exact values. The results indicate we can get reliable values by new method. 

How far is our method useful when the size of the chain is enlarged? Fig. 3(c), 
which shows how ratio P in our method increases for larger A^'s, answers this ques- 
tion. We see that with our present choice of states reliable calculation on A^ > 32 
and n > 4 lattices is difficult when T < 1.0. 

Section 4 Discussions 

In this paper we suggested a new approach for Monte Carlo study of quantum 
spin systems suffering from the negative sign problem. Essential point of our method 
is to choose a set of states for the path integral which is appropriate to numerical 
calculation. Conventional choice of the set, whose states consist of eigenstates of 
z-component of the Pauli matrix on each lattice site, is the simplest one but the 
N.S. problem turns out to be serious with this set. Since any other choice is possible 
as long as the set is complete, we employ a set made of eigenstates of diagonal- 
ized Hamiltonian for every two neighboring sites. We believe this choice is better 
approximation to true ground state and excited states. 

We applied our method to one-dimensional quantum spin 1/2 system with next- 
to-nearest neighbor interactions, which suffers from the N.S. problem in the anti- 
ferromagnetic case. Then remarkable improvements were found as shown in Section 
3. Monte Carlo results also showed, however, that even in our method one would en- 
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counter difficulties for low T on large lattices because the ratio of negatively signed 
configurations increases up to ~ 0.5 when temperature goes down or lattice size is 
enlarged. So we should emphasize that what we present here is not a way to com- 
pletely solve the N.S. problem but a prescription to improve numerical calculations. 

Comments on ergodicity are in order. In quantum Monte Carlo simulations it is 
known that local flips are not sufficient to guarantee ergodic update of the system 
0. One should add some kinds of global spin flips — n-direction global flips to 
change 2;-component of spins, x-direction global flips and transverse global flips to 
change the winding number, and so on — in updating procedure in order to make 
phase space big enough to observe satisfactorily precise results. What global flips are 
necessary and sufficient is not a trivial question, but from technical point of view it 
is preferred to include minimum kinds of global flips which are as simple as possible. 
From our experiences n-direction global ffips seem inevitable for this purpose. As 
results in Section 3 indicate, we can measure system's energy without any other 
flips in wide range of temperature. Transverse global ffips and x-direction global 
flips, however, might be important for low temperatures. It should be investigated 
in future study. 

Although a complete solution for N.S. problem could not be obtained, our results 
show much improvement in numerical study, which inspires us with confidence that 
quantum Monte Carlo methods are effective if a set of states for the path integral 
is appropriately chosen. More improvements would be possible by selecting larger 
cluster of spins to construct states which diagonalize the Hamiltonian. 
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Figure captions 
Fig. 1 

(a) Way to divide Hamiltonian on A^=8 chain in the conventional approach. 
Filled circles, squares, diamonds and triangles indicate terms in Hi, H2, H3 and H^, 
respectively. 

(b) Way to divide Hamiltonian on A^=8 chain in our approach, where odd and 
even sites are renumbered as shown. Open triangles and diamonds represent terms 
in Hi and H2, respectively. Open circles denote terms halved by Hi and H2- 

(c) Examples for twisted n-direction global flips when N — S. Open circles and 
plaquettes with vertical trotter direction denote site locations and elementary four 
spin interactions, respectively. All spins on each path shown by bold solid line should 
be simultaneously flipped. 

Fig. 2 

Energy in the ferromagnetic case versus temperature T. Data are calculated on 
A^=8 chain by diagonalization method (exact) and by the conventional approach for 
trotter number n — 2, 3 and 4. For each T and n twenty thousand sweeps are done 
for the thermalization, one hundred thousand for measurement. 

Fig. 3 

Data obtained by new method in the antiferromagnetic case. Each measurement 
is done using ten thousand conflgurations generated after one thousand sweeps for 
the thermalization. 

(a) Ratio of negatively signed configurations on A^=8 chain as a function of 
temperature. Results from conventional approach are also plotted for comparison. 
In simulations by conventional method twenty thousand sweeps are done for the 
thermalization, one hundred thousand for measurement. 

(b) Energy of the system observed on A'"=8 and n—2, 3 and 4 lattices. 

(c) Ratio of negatively signed conflgurations as a function of the number of lattice 
sites. Trotter number n is kept 4 here. 
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